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We study the real-time evolution of large open quantum spin systems in two spatial dimensions, 
whose dynamics is entirely driven by a dissipative coupling to the environment. We consider different 
dissipative processes and investigate the real-time evolution from an ordered phase of the Heisenberg 
or XY-model towards a disordered phase at late times, disregarding unitary Hamiltonian dynamics. 

The corresponding Kossakowski-Lindblad equation is solved via an efficient cluster algorithm. We 
find that the symmetry of the dissipative process determines the time scales which govern the 
approach towards a new equilibrium phase at late times. Most notably, we find a slow equilibration 
if the dissipative process conserves any of the magnetization Fourier modes. In these cases, the 
dynamics can be interpreted as a diffusion process of the conserved quantity. 

PACS numbers: 03.65.Yz, 05.70.Ln, 75.10.Jm 


I. INTRODUCTION 


Simulating the real-time evolution of large quantum systems is one of the major challenges of modern theoretical 
physics, ranging from condensed matter physics to high-energy physics GH3- Unlike quantum systems in thermal 
equilibrium for which quantum Monte Carlo simulations have proven to be an extremely powerful tool, there is no 
generally applicable approach to far-from-equilibrium quantum systems: Exact diagonalization techniques become 
technically impossible for large systems due to the exponential growth of the Hilbert space with the system size. 
On the other hand, the application of the quantum Monte Carlo method based on importance sampling fails for 
real-time simulations due to a severe sign or complex weight problem. For gapped 1-dimensional systems with small 
entanglement, however, tensor network states underlying the density matrix renormalization group SI provide a 
good basis for simulating the real-time dynamics for moderate time intervals SS3- Moreover, simulations in the 
classical statistical or truncated Wigner approximation are valid in the limit of large occupation numbers which 
appear, for instance, at the early stages of relativistic heavy-ion collision or in the study of non-perturbative fermion 
production [l^ - ll7j . Finally, there have been attempts to apply stochastic quantization to real-time problems based 
on the complex Langevin equation even though the problem of run-away trajectories is not completely settled yet 

"Tit. 


One reason for the complexity of simulating the real-time evolution of large quantum systems on classical computers 
is the fact that isolated quantum systems tend to evolve into complicated entangled states such as those of Schrodinger’s 
cat. For this reason, it was proposed to use specifically designed quantum devices - so-called quantum simulators - 
to mimic quantum systems that are difficult to simulate classically due to the high degree of entanglement [22|. Due 
to the rapid development of atomic, molecular, and optical (AMO) physics in recent years, it has become possible to 
realize quantum simulators in systems of ultracold atoms j23j . As a consequence, applications in atomic and condensed 
matter physics fiulljol as well as in a particle physics context fjoi - fjiil l are widely discussed. 

Furthermore, real quantum systems are usually not isolated but coupled to a dissipative environment resulting 
in decoherence and the suppression of Schrodinger cat states. Hence, it should be simpler to simulate quantum 
systems on a classical computer if we suppress the unitary Hamiltonian evolution and consider a pure, or at least 
sufficiently strong, dissipative coupling to an environment. The corresponding dynamics is then governed by the 
Kossakowski-Lindblad equation [40L |41| , which is the most general non-unitary, Markovian time-evolution equation 
for a density matrix which preserves the basic properties of Hermiticity and positive semi-definiteness. The study of 
the interplay between the unitary Hamiltonian dynamics and the dissipative Markovian dynamics which may result 
in non-equilibrium (quantum) phase transitions has attracted much interest in recent years from both the theoretical 
and the experimental side |42l - l52l . 

In this publication, we investigate the real-time evolution of large open quantum spin systems whose dynamics 
is entirely driven by specific randomized quantum measurements performed on pairs of neighboring spins. The 
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measurement processes give rise to a dissipative coupling to the environment, such that the system is described by 
the Kossakowski-Lindblad equation which can be solved with an efficient loop-cluster algorithm [53j, |54j]. Unlike 
quantum trajectory methods, which are widely used in quantum optics to solve the Kossakowski-Lindblad equation 
and which scale exponentially with the size of the Hilbert space [551 ], we can simulate the real-time evolution of the 
quantum system over arbitrarily long time intervals in any spatial dimension since the computational effort scales only 
linearly with the spatial volume. We then find that the system is driven from an initial equilibrium state to a new 
equilibrium at final times. Compared to a previous study of some of the authors [56 ], we investigate a wider class of 
dissipative processes with different symmetries as well as using a variety of initial density matrices, corresponding to 
the s = 1 Heisenberg model as well as the quantum XY model at low temperatures. We emphasize that the associated 
Hamiltonians are only employed for preparing the initial density matrices but not for driving the subsequent real-time 
dynamics. 

More specifically, we give a detailed account of the equilibration dynamics for the different dissipative processes. 
We show that the symmetry of the specific process determines the time scales which govern the approach towards a 
new equilibrium phase at late times. Most notably, we find a slow equilibration if the dissipative process conserves 
any of the magnetization Fourier modes. In these cases, the dynamics can be interpreted as a diffusion process of the 
conserved quantity. 

This paper is organized as follows: In Secs. Ill Al III Cl we recall the basic equations which govern the real-time 
evolution of a density matrix, most notably the Kossakowski-Lindblad equation, and introduce different measurement 
processes of pairs of neighboring spins. In Sec. IIII Al we discuss the real-time evolution of an anti-ferromagnetic 
Heisenberg model initial state subject to the different measurement processes in great detail and we analytically 
derive expectation values for the new equilibrium state at final times. Conclusions and an outlook are presented in 
Sec. HV] We summarize the rules for forming loop-clusters for the different models and measurement processes in 
Appendix [A] The real-time dynamics of a ferromagnetic Heisenberg model initial state and of a quantum XY-nrodel 
initial state are discussed in Appendix [5] and Appendix [Cj respectively. 


II. MEASUREMENT-DRIVEN DISSIPATIVE DYNAMICS 

The real-time dynamics of a closed quantum system from an initial time to to a final time tf is governed by the 
unitary time-evolution operator 


U(tf,t 0 ) =exp(-iH(t f -t 0 )) = U^(t 0 ,t f ) , (1) 

where we assumed a time-independent Hamiltonian operator H. Starting from an initial state which is specified by a 
(not necessarily thermal) density matrix p(to), the expectation value of any observable O at a later time is given by 

(0)(t) = Tr [p(t)0] = Tr[p(t 0 )U(to,t)OU(t,t 0 )\ , (2) 

where we employed the von-Neumann equation 

p{t) = U(t,t 0 )p{to)U(t 0 ,t) . (3) 

The expression ([2]) may be rewritten as a real-time path integral along the Schwinger-Keldysh contour [57l [58| or, 
assuming a thermal density matrix p(to), along the Konstantinov-Perel’ contour [59[ . The fact that the matrix elements 
of the time-evolution operator are in general complex leads to a severe sign problem, preventing the application of 
quantum Monte Carlo techniques. Accordingly, we will consider a simpler situation in the following. 


A. Real-time evolution by sporadic measurements 

We assume that certain observables Ok are measured at intermediate times tk with k £ {1,..., V}, resulting in 
measurement results Ok■ The individual measurement results are associated with projection operators P Qfc , projecting 
on the subspace of the Hilbert space which is spanned by the eigenvectors of Ok with eigenvalue Ok■ The projection 
operators are Hermitean P 0k = idempotent P% k = P 0k , and fulfill J2o k P°k = ^ We note that, given a certain 
density matrix p(tk) at time tk, a measurement of the observable Ok with result Ok yields a new density matrix 

jn \ PoMPo. 

P l tk) = - 7 —r- 


( 4 ) 
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where p(o k ) = tr [p(t k )P 0k \ denotes the probability of obtaining the measurement result o k . Alternating unitary real¬ 
time evolutions according to von-Neumann’s equation ([3]) and quantum measurements according to ((H) for a given 
initial density matrix p(to) yields after N steps 

p(ijv) = —7 - rP ON U(tN,tN-l) ■ ■ ■ Po 2 U(t 2 ,ti)Po 1 U(ti,to)p{to)U(to 1 ti)P 0 - l U{ti 1 t 2 )Po 2 ' ' ' tjv)-P OJ v • 

( 5 ) 

Here, p(oi, 02 , ■ ■ ■, on) denotes the probability of finding the measurement results {cq, 02 ,..., On} upon starting with 
the initial density matrix 

P(t 0 ) = , (6) 

i 


with 0 < pi < 1 and ^jPi = 1. Accordingly, the probability of reaching a certain final stat e / ) at the final time tf 
after a sequence of N measurements with measurement results {01,02,..., ojv }, is given by [60( 

Pp 0 ,f(.Oi,O2,... 1 O N ) = (f\U{tf,t N )p(t N )U{tN,tf)\f) = 

y^ j Pi(i\U(to,t 1 )P 0l U(t 1 ,t 2 )Po 2 ■ ■ ■ P ON U{t N ,tf)\f)(f\U(tf,t N )P ON ■ ■ ■ P 02 U{t 2l ti)P 0l U{t\,t 0 )\i) . (7) 

i 

This transition probability is still plagued by a severe sign problem. In order to make the problem amenable to 
Monte Carlo importance sampling, we further assume in this publication that the real-time dynamics of the quantum 
system is entirely driven by the measurement process, i. e. U(tk+i, tk) = 1. The more challenging step of the real-time 
evolution being driven by both a non-trivial time-evolution operator and by quantum measurements is currently under 
investigation but beyond the scope of this publication. It has to be emphasized, however, that even this simplified 
problem results in interesting highly non-trivial quantum dynamics (5(|, as discussed in more detail in the following 
sections. 

In this case, the probability of reaching a final state | /) at the final time tf = tjv after a sequence of N measurements 
with measurement results {01,..., on } is calculated according to 

P P0 j{oi,...,o N ) =Y^ / Pi( i \ P oi ■■■Po N \f)(f\Po N ---Po 1 \i) ■ ( 8 ) 

i 

In order to derive a real-time path integral along the Schwinger-Keldysh contour leading from the initial time to to 
the final time tN and back, we insert complete sets of states between the individual projection operators. Denoting 
these complete sets by {|nj.)} on the forward branch and by {|ufe)} on the backward branch, with k £ {1, ...,1V— 1}, 
the resolution of the identity operator reads Y^n k \ n k)(rik\ = Yin' | n fe)( n fe| = Accordingly, we find 

N 

Pp 0 ,fi.Oi,...,O N ) = 5>£- £ n {n k -in' k _ l \P 0k <g> P* k | n k n' k ) , (9) 

i n\riN — i,n' N1 k—1 

where we introduced the notation 

(n k - 1 n k _ 1 \P 0k ®P* k \n k ri k ) = (n k -i\P 0k \nk){n' k ^\P 0k \n k )* , (10) 

with (u-oUqI = (ii\ and \nNn' N ) = \//). We note that the doubled Hilbert space of states {\n k n' k )} encompasses 
both branches of the Schwinger-Keldysh contour. 

In the following, we will be mainly interested in the total probability of reaching the final state |/) irrespective of 
the intermediate measurement results {01,02,..., on }, which is calculated according to 


N 


Ppo,f =£---£Ppo,/(°i>---> 0 iv) £ "■ £ n ( n k-in' k _ 1 \P k \n k n' k ) . (11) 

i n\ ,n^ tin _ i,n' N _ 1 k= 1 


01 On 


Here, 


^ = £ P, k ® pi 


( 12 ) 
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is obtained by summing over all possible measurement results Ok at time tk- The final density matrix p(t]\r), which is 
the result of the real-time dynamics generated by the sequence of N measurements, is then given by 

p(*iv) = E-E Po N p 01 p(to)p 01 ■ ■ ■ Pos > (13) 

Oi On 

with Tr [p(t N )] = T,fPpoJ = !■ 

B. Lindblad evolution 

In the previous section we described the real-time evolution of the density matrix if the dynamics is entirely driven 
by a sporadic measurement process at discrete times tk- In the following, we consider the continuous time limit 
tk+i — ffc = e —> 0 and suppose that the measurement happens only with a certain probability per unit of time. 
Accordingly, the quantum system can then be considered as being continuously monitored by the environment. 

This situation can be described by the Kossakowski-Lindblad equation [4CJ, |41[ , which is the most general non¬ 
unitary, Markovian time evolution equation of a density matrix which preserves the basic properties of Hermiticity 
and positive semi-definiteness. It is characterized by a set of Lindblad operators, describing the possible quantum 
jumps the system may undergo at any instant of time 

L o k = VnPo k , ( 14 ) 


which obey 


(l-ey)! + J2 L l k L 0k = 1 . 

Ok 


(15) 


The parameter 7 determines the probability per unit time that a certain quantum jump occurs. Using the notation 
of the previous section, a certain quantum jump operator L 0k corresponds to the measurement of an observable Ok 
yielding the measurement result Ok - 1 The Kossakowski-Lindblad equation is then given by 




L o k p(t)Ll k 


2 L l k L o k p{t ) 



= 7 J2[Po k P(t)Po k 

Ok 


P(t)] 


(16) 


We emphasize the absence of a unitary time evolution term —i[H,p(t)] on the right-hand side as we assume that the 
real-time dynamics is entirely driven by the dissipative process. 

We translate this continuous-time equation back to discrete times since we employ a discrete-time algorithm in our 
numerical simulations. We note, however, that the systems could also be simulated directly in continuous time & 
To this end, we recursively perform a forward finite-difference discretization of the time derivative to obtain 


P(tk) = (1 - e7)p(t fc _i) + £7 E P o k p(tk-i)Po k , 

Ok 


(17) 


where we disregard higher discretization errors. Here, k € {1, ...,7V} labels the discrete times tk at which the quantum 
system potentially interacts with the environment. The recursive equation QU) allows us to express the final density 
matrix /o(tjv) as a function of the initial density matrix p(to) according to 

N 

p(t N ) = (1 -£7) jv p(f 0 ) -(-^(eq)^ 1 - £7)^"* E E " ‘ E P °u '' ’ P °n^o)-P 0il • • • Po (a ■ (18) 

s=l Oi s 

The individual terms have simple interpretations: The first term ~ (1 — e^) N correspond to a time history without 
any interactions whereas the remaining terms all include measurements. In fact, the order s determines how often the 


1 In general, the quantum system may interact in different ways with the environment so that we have to include several different 
observables Ok at any instant of time. For simplicity, we restrict ourselves to a single observable in this section and discuss the more 
general situation in Sec. mm 
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quantum system interacts with the environment. Most notably, the last term ~ (cj) N corresponds to a time history 
with quantum jumps at every discrete time step. 

Following the same construction as in the previous section, we can again derive a real-time path integral along the 
Schwinger-Keldysh contour leading from the initial time to to the final time t n and back. Most notably, the total 
probability of reaching a final state |/) irrespective of any possible intermediate measurement results is given by 

N 

Ppo,f = Yl Pi n (nk-xn' k _ x \{l-e^)t®t + e^P k \n k n' k ) . (19) 

i 711,71^ njv_ k—1 

As a matter of fact, we recover the result for the real-time evolution due to sporadic measurements (|11[) in the limit 
e 7 =l, corresponding to a measurement at every time step with probability 1. 


C. Measurement processes 

In this study, we investigate 2-dimensional systems of quantum spins s = In order to generate dissipative 
dynamics, we apply a given measurement process corresponding to a specific observable Ok- As we will see, different 
measurement processes result in different dynamic behavior. 

In order to set the stage, we start with the simple situation of just two quantum spins at positions x and y. At 
the end of this section, we will generalize the two-spin system to a large 2 -dimensional system, appreciating that 
this generalization is feasible for any dimension. In the following, we take the 3-direction as the quantization axis 
and denote the two eigenvalues of S x by s x = t = \ and s x = = —1. The Hilbert space of the two-spin system is 
composed of four states and denoted by 

^ = {|sx^)} = {|tt),|tl),Ut),|U)} ■ (20) 

1. Measurement process: S 2 

The first measurement process under consideration corresponds to the total spin of the system 

0 (1) =s 2 = (S x + Sy) 2 = S% + S 2 + 2 s x ■ Sy = + 2 S x ■ Sy . ( 21 ) 

This operator conserves each of the components of the spin vector S = S x + S y so that 

[S 2 ,S i ] = 0. ( 22 ) 

The simultaneous eigenstates of S 2 and, for instance, its 3-component S 3 = S x + S 3 are denoted by |S'5' 3 ). They are 
given by the singlet state 


00 > = 7f (I ) - | It » 


(23) 


and the triplet states 


{|11> = | tt > , |lO) = ^(|U) + Ut» , |l-l) = |u)} . (24) 


The projection operators on the different measurement results S = 0 and S = 1, respectively, are then given by 


P 0 = |00)<00| , 

Pl = | 11 )( 11 | + | 10 )( 10 | + |1 - 1)(1 - 1 

In the basis (12011 . these projection operators have the matrix representation 



(° 

0 

0 

°\ 

/2 0 

0 

°\ 

1 

0 

1 

-1 

o ft 1 

0 1 

1 

0 I 

2 

0 

-1 

1 

0 ’ Pl 2 

0 1 

1 

0 


\0 

0 

0 

0 / 

\0 0 

0 

2 / 


(25a) 

(25b) 


( 26 ) 
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We note that certain matrix elements for measuring S' = 0 in the doubled Hilbert space of states | n^n^.) = 
| s x,kS y ,ks' x k s' y k ) are negative 


^x,k^y,k^x,k^y,k | -^0 ® ^0 |s a;> /c+l'Sy,fc-|-l'Sx,fc-t-:iAi/,fc+l) 

1 


^ ^ I..S™ l. i i 5.,.. l, i i 


+ ’ S </,fe + l Us x,fc’ S i,fc + l^ S H,fe’ S a,fc + l ^ Sx , k ’ S y, k +1 ^s s ,fc> s *.fe + l ^ s L,fe’ S B,fc + l ^ S B,fe> s i,fc + i 


-<5. 


-tf. 


*x,fc>Sx,fc+l -Si/.fe s xk ,s yk+1 [J s yk ,s xk + 1 a Xtk ,Sy :k+ 1 s B|fc ,s Xifc+1 > s Xi fc+l s H ,fc ■ ■ s „,fc-t-i 1 ’ 


.) 


(27) 


giving rise to a sign problem in the corresponding real-time path integral © ■ On the other hand, the matrix elements 
for measuring S = 1 are always non-negative 

( s x,k s y,k s x,k s y,k\Pl ® ^1 | s x,fc+l s J/,fc-|-l s x,fc+l S 2 /,fc+l) = 

4 ^ ^ S *,fc' S *,'= + l^ S !/,'=’ S »,fc + l^ s L,fc> S x,fc + l < ^ S l / ,fc’ S (,lc + l ^ Sx < k ’ S y- k + 1 ^ S y- k ’ Sx < k + 1 ^ S 'a,,k’ S y,k+l^ S 'y,k’ S 'a 1 ,k + l 

^ r ^ Sx ,k^ :Cj k+l^ay : k,ay : k+l^S x k ,s' y k+1 ^s' y k ,S x k+1 T ^S xk ,S H ,fc + l ^ s y,k ,a x ,k+l ^ s x , k > s L,fc + l ^ S 'y,k < S y,k + l ) ' 

Remarkably, the sign problem arising from m is eliminated by averaging over both measurement results 


( s x,kSy t kS X 'i~Syj,\Po <S> P 0 + -Pi <8> P\ | s x,k+lSy,k+lS x ^k+\ s y ,fe+i) — 

3( 4 


s x, k t s x, k +l^ay ik ,Sy ik+1 ^s xk ,s xk+1 ^s yk ,s' yk+1 T 3s X ' k ,Sy ik+ i ^S yk ,s X ' k+1 ^s xk ,s yk+1 ^a yk ,s' xk+1 ) — 0 > 


(29) 


so that the corresponding real-time path integral CD is not plagued by a sign problem anymore. 


2. Measurement process: S x Sl 

The second measurement process under consideration corresponds to a measurement of the products of 1- 
components of the quantum spins 

0 (2) = SlSl . (30) 

Accordingly, this measurement distinguishes whether the 1-components of the two spins are the same or different. 
Equivalently, we could have also chosen a measurement process S x Sy regarding the 2-components of the two spins. 
As we choose the 3-direction as the quantization axis ©, the two eigenstates of O^ corresponding to parallel spins 
are denoted by 


{l II 1 > = 75 (I tt ) + | 44 » , |ll2> = ^(|U> + Ut»} » ( 3 !) 

whereas the two eigenstates corresponding to anti-parallel spins are given by 

{Ui) = 7f(|tt)-|U» , I-H2 > = 75 (| tl) - Ut»} • ( 32 ) 


The projection operators on parallel and anti-parallel 1-components of the two spins, respectively, are then calculated 
according to 


iii H Hi X Hi I +1II 2 X h\ 

Px = \h)(h\ + \h)(h\ 

Again, expressing these projection operators in the basis (12U1) . we obtain 



/l 0 0 1\ 


/ 1 0 0 -1\ 

1 

0 110 

5 1 

0 1-101 

2 

0 110 

’ P#= 2 

0-11 0 


\1 0 0 1/ 


\-l 0 0 1 / 


(33a) 

(33b) 


( 34 ) 
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Like before, certain components of ar 
On the other hand, the matrix elements 
have 

k s y,k\Py; ® Pft\ Sx , k +l s y 


^ ^ ^ s x,k > s x,fc + l ^ s y,k 5, fc-t-1 ^ s x t k 


s y,k+l s x ,k+l S y,k+l 


L> = 


( s x,kSy,kS x ,k s y , 

1 / 

^ S 'x,k’ S x,k+l^ S y,k’ S 'y,k+l + S *x,k-*x,k+ l^Sy 




,5 y,fc + l ^ 






<5,' < 


(Sx,fcSy,fc'S a , ife S yife |P|| <S> Pj| |Sx,fc-K s 2 /,fc+l S x,fc+l S y,fc-|-l) : 

1 ^Sy,k ,Sy,fe+l ^ S x, k ’ S 'x, k + 1 ^ S y,k ’ S y,k + 1 ,k s ®, fc+1 ^ s y, k , ~ Sy, k + 1 ^ s ' x ,k S x ,k + l ^ S y,k’~ S y,k + 1 


^ ^ ^ s x,fc, s x,fc + 


+ ^,fe,Sx,fe+l^Sy,fe,Sy,fe+l^< ife ,-s' X|fe+1 ^ ifc -S^ fe+1 


(35a) 


(35b) 


As before, the sign problem arising from (135al) is completely eliminated by averaging over both measurement results 
(sx,kS y ,ks' Xik s' ytk \P\\ ® Pf + P#® P^\sx,k+is v ,k+is' Xtk+1 s' y k+l ) = 

2 ( ^ Sx - k ’ Sx :k+l^y,k,Sy : k+l^s' xk ,S l xk+1 ^s l yk ,s l yk + 1 T ^Sx,fe , — **,4+1 ^ s V,k Sy,k+1 ^ s ' x ,k s i,k + l ^ S 'y,k ’ ~ S y,k+1 ) — 9 ’ 


3 . Measurement process: S x S + + S x S y 


Finally, the last measurement process under consideration corresponds to the observable 

o (3) = S+S+ + S-S- = 2 - S 2 S 2 ) , (37) 

where we introduced the raising and lowering operators according to S x = Sj. ± iS 2 . The intriguing feature of this 
observable is that it conserves the difference of the 3-components of the two spins 

[S+S+ + S-S~,S 3 X - P 3 ] = 0 . (38) 

This process measures the correlation between parallelism and anti-parallelism of the 1-components and 2-components 
of the two quantum spins. Given that the 1-components are the same whereas the 2-components are different, the 
eigenvalue of O^ is +1 and the corresponding eigenvector is given by 

I + ) = 72 (I ^ ) + I ^ )) • ( 39 ) 

In the opposite case, where the 1-components are different but the 2-components are the same, the eigenvalue of O^ 
is —1 and its eigenvector reads 

|-) = 7f(|tt)-|U» . (40) 

Finally, there are two degenerate eigenvalues 0 corresponding to the case where the 1-components and 2-components 
are both either the same or different. Because of (1551) . we can further distinguish these states by their simultaneous 
eigenvalue of S 3 — S 3 and we denote them by 


(|o + > = |U>, |o_) = Ut)} • 


The three projection operators on the different measurement results are thus given by 

P + = I + )( + I > 

P 0 = |0+)(0+|+|0_)(0_| , 

p - = l-X-l • 

Equivalently, the matrix representation in the basis reads 



/I 0 0 1\ 


/0 0 0 0\ 

1 

0 0 0 0 

5 

0 10 0 

2 

0 0 0 0 

) Po — 

0 0 10 


\1 0 0 1/ 


Vo o o o/ 



/ 1 0 0 
0 0 0 
0 0 0 
\-l 0 0 



(41) 

(42a) 

(42b) 

(42c) 


( 43 ) 







Like before, the negative entries in P_ give rise to a sign problem in the real-time path integral 
remaining matrix elements are non-negative 

( s x,kSy,kS x ,k s y,k\P- <8> P-\s x ,k+lSy t k+lS x j i . + 1 S y<k+ i') = 

/ ^ s x,k s x,k+is xk s x k+1 5 SxyctSy j c S Sxk+ltSy j c+1 5 s ^ k ^ k S s ^ k+itS ^ k+i , 

( s x,k s y,k s x,k s y,k\P+ ® P+\ s x,k+lSy,k+lS Xt k+l s y } k+l) = 

^S Xlk ,Sy,k d Sx k + 1 ,Sy tk + 1 ^s x k ,S y k ^s’ x k+ 1 ,s' y k+1 ) 

{s x ,kSy^kS x k S y k \Po (8) Pq |s a ; ) A,+iSy i fc+iS a , jfc+1 S yifc+1 ) = 

(1 + 4s a , i fcS a , t fc+l)(l + 4:S Xtk s' x j, +1 ) x A A 

4 Sx,k,-Sy,k Sx,k+ 1: -Sy,k+1 s ' X 'k’- S y,k S 'x,k+1 > ~ S 'y,k+l 

The sign problem arising from (144al) is again completely eliminated by averaging over all three measurement results 
( K s x,kSy,kS xk S yk \Pj r (gl P + + Pq ® Pq + P_ ® P_ | s x,k+lSy,k+\S xk ^S yk+ ^ ) = 


whereas all 

(44a) 

(44b) 

(44c) 


^ s x,kS x ,k ( s x,kS x k + S Xj / c +iS. J , jfc _|_ 1 ) ^s Xtk ,s y ,k^s x ,k+r>Sy,k+i^s' xk ,s' yk ^s' xk+1 ,s' yk+1 
4Sa;,fcS Xl fe(Sa:,A: + S x ,fc-|-l) (s x k + S X] fc + i) ^s Xtk ,-s ytk Sa Xt 


k + l, s y,k + l^ s Xik , + S (,fe + 1 — ^ 


(45) 


4- Generalization to large higher dimensional systems 

We now generalize the results for the two-spin system to large 2-dimensional systems of quantum spins s = 4 
on bipartite square lattices of the size L x L with periodic boundary conditions. In fact, the generalization to any 
number of dimensions is feasible and straightforward. To this end, we assume that any measurement process affects 
only nearest-neighbor quantum spin pairs. Accordingly, we can take advantage of the results of the previous section 
for the two-spin system. 

Specifically, in order to allow for an efficient implementation of the real-time evolution, we discretize the system in 
the following way (cf. Fig. |TJ): In a first step, all pairs of neighboring spins separated in the 1-direction at x = ( Xi,X 2 ) 
and y = (x\ +1, X 2 ) with even x\ can be measured. The second step allows for measuring spin pairs in the 2-direction 
at x = ( 21,22 ) and y = ( X\,X2 + 1) with even X2 ■ In the third and fourth step, the measurement process affects spins 
with odd X\ and odd X2, respectively. Accordingly, all possible interactions between nearest-neighbor spins occur 
during these four steps, which can then be repeated an arbitrary number of times. 

For a sporadic measurement process, all of the neighboring spins are measured simultaneously according to the 
four-step scheme. The corresponding probability of reaching the final state |/) irrespective of the intermediate 
measurement results m is then given by 

N 

PpoJ = 5>£- £ IIII < s " -l s y,k-lS xk _iS y k _ 1 \Pk,xy\s x ,kSy y kS x k s y t k) > (46) 

i ni.n'j n N - 1 ,n' N _ 1 k=l {xy) 

where n k = [s x ,fe] represents the 2 N2 possible spin states. We emphasize that the 27V 2 projectors P k , xy are chosen 
according to the aforementioned four-step scheme, i. e. for k € {1, 5, 9,...} only those 7V 2 /2 projection operators are 
used which connect neighboring spins separated in the 1-direction with even x±, and likewise for the other three steps. 

In the following, we study the more natural situation where the measurement process does not affect all quantum 
spins simultaneously but only a small number of randomly chosen neighboring spins. This situation is again described 
by the Kossakowski-Lindblad equation, which now results in 

N 

ppoj = ££ £ ••• £ n n <^- i ^, fc - i 4, fc - 1 4, fc - 1 |( 1 -7) 1 ® 1 + e7 p fc , x ,| Sx , fc ^4, fe 4, fe ) • ( 47 ) 

i ni,nj n N - 1 ,n ' N _ 1 fc= 1 < xy) 

In fact, we retain the aforementioned four-step scheme for discretizing the system for our numerical simulations. We 
note, however, that four discrete time steps are considered as one physical time step in the spirit of a Suzuki-Trotter 
decomposition such that t k ^ = e^k/A. Moreover, we emphasize that the particular order of the four-step scheme is 
irrelevant for the Lindblad evolution since the interaction is anyway randomized in the continuous time limit. By 
taking the limit ey —> 0, only a small number of randomly chosen neighboring spin pairs interact at any instant of 
time. 
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(1) 


( 2 ) 


£2 



Xi 


FIG. 1: [Color online] Visualization of the discretization of the physical system. The oval-shaped objects denote the nearest- 
neighbor spin pairs which may interact at discretization step 1-4, respectively. 


The physical process under consideration is the following: Starting from a thermal initial density matrix p(t $), the 
measurement process induces a non-trivial real-time evolution, resulting in a final density matrix p(tAr). In order to 
define p(fo), we will choose different model Hamiltonians H corresponding to the s = i anti-ferromagnetic (AFM) and 
ferromagnetic (FM) Heisenberg model as well as the quantum XY-model. The Euclidean time interval [0,/3], where 
/? = 1/T is the inverse temperature, together with the real-time interval [0,tjv] then forms a closed Konstantinov- 
Perel’ contour in the complex time plane. We again emphasize that H is only used to prepare an ensemble of initial 
states whereas the real-time evolution is entirely driven by the dissipative measurement process. In order to simulate 
this process, we will employ a multi-cluster algorithm [531 . |5j|. We summarize the rules for forming loop-clusters for 
the different models and measurement processes in Appendix [A] 


III. RESULTS AND DISCUSSION 

We now come to the results regarding the dissipative measurement-driven dynamics of spin models, which are 
based on numerical simulations employing a loop-cluster algorithm. We study in detail the dynamics for different 
measurement processes as summarized in Tableland show the crucial role of the symmetry of the measurement process 
for the real-time evolution. In the following, we investigate the initial state corresponding to the antiferromagnetic 
Heisenberg model in great detail. The initial states corresponding to the ferromagnetic Heisenberg model as well as 
to the quantum XY-model are discussed in Appendix [Bl and Appendix [Cl respectively. 


A. Heisenberg anti-ferromagnet initial state 

In this section, we consider an initial density matrix po corresponding to the anti-ferromagnetic Heisenberg model 
<E3]h where we choose the 3-direction as the quantization axis. The ensemble of initial states is then prepared by 
means of the Euclidean-time cluster rules (|A3jl . In the following, the 3-component of the staggered magnetization 






measured observable 

S 2 = (S X + Sy) 2 

si Si or sis' 2 

C »1 nl q 2 o 2 

conserved Fourier mode S(p) 

(0,0) 


( 7 T, 7 r) 

final equilibrium values A(p) 

Eq. (1611) 

Eq. ([63]) 

Eq. E]) 


TABLE I: The measurement processes under consideration, along with their conserved Fourier mode and their final equilibrium 
values A(p). 
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order parameter is denoted by 


M s = ^2(-l) Xl+X2 S 3 x . (48) 

X 

In order to analyze the subsequent real-time dynamics, it is useful to introduce the Fourier modes for the two- 
dimensional square lattice of size L x L according to 

S(p) = ^ exp (ipx) ^ exp {ip\X\ + W 2 X 2 ) Sf r . (49) 

X X 

Accordingly, the 3-component of the staggered magnetization will be denoted as the ( 7 r, 7r)-mode whereas the 3- 
component of the uniform magnetization 


M = J2 S * 

X 


(50) 


is denoted as the (0,0)-mode. 


1. Dynamics for different measurement processes 

Starting with an initial ensemble at low temperature /3J = 5L/2a = 40 on the square lattice of size L x L with 
L = 16a, where a is the lattice spacing, we investigate the dissipation-driven dynamics for different measurement 
processes. We emphasize that the various measurement processes strongly differ in their conservation properties. On 
the one hand, the measurement process O W conserves the (0,0)-nrode whereas the measurement process O^ results 
in the conservation of the (zr, 7r)-mode. On the other hand, the measurement process O^ does not conserve any 
Fourier component. 

Considering the Lindblad evolution with eq = 0.05 for different measurement processes, we display the real-time 
dynamics of a variety of Fourier modes 

(\S(p)\ 2 )(t k ) = Tr[p(t k )\S(p)\ 2 ] , (51) 

with k £ {0,..., IV} in Figs. [2] 0] We checked the dependence on varying ej to guarantee that we are effectively 
simulating a continuous Lindblad process with the discrete-time algorithm. Most notably, the different conservation 
properties of the measurement processes are reflected in the time-dependence of the Fourier modes. While the 
conserved quantities i.e. the (0, 0)-mode for O W and the ( 77 ,7r)-nrode for - do not equilibrate at all, the Fourier 
modes in the vicinity of the conserved quantity show a much slower equilibration rate than Fourier modes remote 
from the conserved quantity. On the other hand, all Fourier modes show a rapid equilibration for the measurement 
process OO) for which there are no conserved Fourier modes present. 

After an initial phase, which is studied in more detail in the next section, the various Fourier modes approach their 
ultimate new equilibrium exponentially 

(IS( P )\ 2 )(t) ‘ ~°° A( P ) + B(p) exp (-i/r(p)) , (52) 

where 1/V(p) denotes the equilibration rate and A(p) is the final equilibrium value. Employing the definition (1491) . we 
find 


(I‘S'(P)| 2 )W =^2^2exp{ip(x-y))(S^)(t) , 

x y 


(53) 


so that 


/V 4 

^(|^(p)| 2 )(t) = iV 2 ^((^) 2 )(t) = x , (54) 

p X 

where N denotes the number of lattice points in each spatial direction such that L = Na. Accordingly, in order to find 
A(p), we need to calculate the spin-spin correlation function (S^Sy) in the final state. In fact, the final density matrix 
p(t —> 00 ), to which the system is driven by the measurement process, is constrained by the conserved quantities of 
the measurement process and is proportional to the unit matrix in each symmetry sector. 
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FIG. 2: [Color online] (Heisenberg anti-ferromagnet initial state) Time evolution of certain Fourier modes (|<S'(p)| 2 )(t) for the 
measurement process CAR Left: Linear plot for a short time interval, cf. also Fig. lc in j5§|. Right: Log-log plot for a long 
time interval. The error bars are of the order of the symbol sizes and the lines are included to guide the eye. We initialize an 
anti-ferromagnetic Heisenberg model at low temperatures such that the (7r, 7r)-mode is large whereas the (0,0)-mode vanishes. 
The remaining parameters are 47V T = 512, L = 16a, /3J = 5L/2a = 40 and ey = 0.05. In order to determine the time evolution, 
we performed 10 6 Monte Carlo measurements. The horizontal line corresponds to the analytically derived final equilibrium 
value (l6lbl) . We emphasize that the (0, 0)-mode is exactly conserved during the time evolution. 


In Fig. [2l we show the real-time evolution of a variety of Fourier modes for the measurement process . As 
already discussed, we find that the (0, 0)-mode is exactly conserved during the time evolution. On the other hand, 
the Fourier mode which is closest to the (0,0)-mode has the slowest equilibration rate whereas the Fourier modes 
which are further away equilibrate more rapidly. As a matter of fact, all Fourier modes except the (0,0)-mode are 
then driven towards the same equilibrium value AA>(p ^ (0,0)), which will be discussed below in more detail. 

The detailed equilibration, however, is rather intricate since a non-trivial attractor A(t) is formed which approaches 
A^\p ^ (0,0)) for t —> oo: Due to the fact that the slowest Fourier mode approaches A^^(p ^ (0,0)) from below, 
all other modes finally reach it from above due to (1541) . This, however, has far reaching consequences for the time 
evolution of the single Fourier modes. On the one hand, Fourier modes which start off above the attractor A(t) fall 
onto it from above so that their amplitudes decrease monotonically in time. On the other hand, Fourier modes which 
start off below the attractor A(t) are enhanced in a first step so that they are driven towards the attractor from 
below. Once these Fourier modes reach the attractor, their amplitude again evolves along the attractor and decreases 
monotonically in time. As a consequence, all Fourier modes except the slowest one will at some point in time fall onto 
the attractor, which always lies above A (1 \p ^ (0,0)). 

Moreover, it turns out that the time needed to approach the attractor A(t) in the first instance depends on the 
specific Fourier mode. As time evolves, more and more Fourier modes approach the attractor and, subsequently, show 
the same relaxation dynamics. In the end, all Fourier modes except the slowest one lie on the attractor such that the 
final approach towards the new equilibrium at t —> oo is determined by the equilibration rate l/r(p) of the slowest 
Fourier mode. 

In order to derive an analytic expression for A (1 \p ^ (0,0)), we consider the system to be prepared at low 
temperature T — > 0 at the initial time such that 

A (1) (p = (0,0)) = (| S(p = (0,0))| 2 )(f o ) ^ 0 , (55) 

corresponding to the ground state of the system for which the uniform magnetization vanishes. Due to the fact that 
the (0, 0)-mode is exactly conserved during the time evolution, the system is driven to a final equilibrium ensemble 
for which the 3-component of the uniform magnetization still vanishes. In general, the final state partition function 
in any sector of magnetization M £ {—TV 2 /2, ...,7V 2 /2} is calculated according to 

2 TT 

Z M \J] = U exp (^2^ S z) S J: z sim = ^ f dX Y{ ex P (^Z( iX + jz) S z- iXM ) ■ (56) 

X 03 = 4-i 2: / X C3 = 4-I 2 

Restricting ourselves to the M = 0 sector of the Hilbert space according to the chosen initial state, we find that there 
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FIG. 3: [Color online] (Heisenberg anti-ferromagnet initial state) Time evolution of certain Fourier modes (|S'(p)| 2 )(t) for the 
measurement process O ^ on a linear plot. The error bars are again of the order of the symbol sizes and the lines are included to 
guide the eye. The parameters are as in Fig. [2] and the horizontal line corresponds to the analytically derived final equilibrium 
value (EH) . We emphasize that all Fourier modes equilibrate very rapidly. 


are 


2tt 


Z M = o[0] = d- f d\ [2 cos (|)] N = 


N 2 '. 

N 2 \N 2 ] 
2 • 2 - 


(57) 


states which are equally probable. The spin-spin-correlation function is then calculated according to 

2tt 

(S 3 X S 3 V ) = 1 


2ttZm=o [ 0 ] 


dA II E S *Sy™ P(E 1 ^)- 


(58) 


o * s 3=±| 

We trivially obtain i for x = y. On the other hand, for i/j/we find 

2 7T 


<S ' S » 3) = (l) = - 


4(fV 2 -l) ‘ 


Accordingly, the spin-spin-correlation function can be written as 

N 2 5 x ,y - 1 


(S d x S. } = 


4(iV 2 -l) ‘ 

As a consequence, upon evaluating (l53l) we find that the Fourier modes take the final values 

A«(p=(0,0)) = 0, 

N 4 


AW(p?(0,0)) = 


4:(N 2 - 1) ‘ 


(59) 

(60) 

(61a) 

(61b) 


The measurement process O^, on the other hand, has the simplest final state as it does not conserve any Fourier 
mode. As a consequence, all Fourier modes equilibrate very quickly as shown in Fig. [3] We note that the equilibration 
of the distinct Fourier modes is not completely independent from each other because of (1541) . which gives rise to an 
overshooting of several Fourier modes, e. g. the (57r/8, 57r/8)-mode in the figure. However, we do not observe the 
occurrence of a non-trivial attractor as there is no slow mode present in the system. In fact, the final equilibrium 
corresponds to a state where each of the 2 N spin states is equally probable. Therefore, the spin-spin-correlation 
function in the final equilibrium is easily calculated according to 




* S 2 =±i 


o3 c3 _ _ r 
°x°y ~ J 


(62) 
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FIG. 4: [Color online] (Heisenberg anti-ferromagnet initial state) Time evolution of certain Fourier modes (|<5'(p)| 2 )(t) for the 
measurement process 0 (3 b Left: Linear plot for a short time interval. Right: Log-log plot for a long time interval. The error 
bars are again of the order of the symbol sizes and the lines are included to guide the eye. The parameters are as in Fig. [5] and 
the horizontal line corresponds to the analytically derived final equilibrium value (tTlbl) . We emphasize that the (7r, 7r)-mode is 
exactly conserved during the time evolution. 


so that all Fourier inodes take the same final value 


A m {p) = — ■ ( 63 ) 

Finally, we display the real-time evolution of the Fourier modes for the measurement process O*- 3 - 1 in Fig. 0] As 
already discussed, this measurement process conserves the ( 7 r, 7r)-mode, corresponding to the order parameter of the 
system which is large at low temperatures T — > 0. In fact, finite-volume chiral perturbation theory predicts 

Ad 2 ! 4 / r \ n 

A^ ) {p= (7T,7T)) = (\S(p= (7T,7r)| 2 )(f 0 ) ~- ^—J2 an yyLj ’ ( 64 ) 

where the first three coefficients are given by ao = 1, ai = 0.45157 and 02 = 0.082803 [6^464| . The low-energy 
parameters are the staggered magnetization density M s = 0.30743(l)/a 2 , the spin stiffness p s = 0.1808(4) J and 
the spin-wave velocity c = 1.6585(10) Ja [fy|, (66|. Accordingly, the typical length scale is given by £ = c/(27rp s ) = 
1.459(3)a. 

Due to the fact that both measurement processes O W and O^ have a conserved mode, the real-time evolution of 
the Fourier modes shows several similarities, most notably the occurrence of a non-trivial attractor A(t). Nevertheless, 
there is also an important difference: The Fourier mode with the slowest equilibration rate towards the final equilibrium 
value A^ 3 \p ^ ( 71 , 71 )) is the one which lies closest to the ( 71 ,7r)-mode. As a consequence, the slowest Fourier mode 
approaches A^ 3 \p ^ ( 71 ,7r)) from above whereas all other modes finally reach it from below due to (1541) . Accordingly, 
the corresponding attractor A(t) lies below A^ 3 \p ( 71 , 71 )) and increases monotonically in time. It has to be 

emphasized, however, that the attractor’s distance from the final equilibrium value A^ 3 \p ^ ( 71 , 71 )) at early times is 
much larger than for the measurement process due to the larger amplitude of the slowest mode. In similarity 
to the measurement process however, all Fourier modes except the slowest one will, at some point in time, fall 
onto the attractor such that their final approach towards equilibrium is determined by the equilibration rate of the 
slowest mode. 

In order to derive a simple approximation for the final equilibrium value A^ 3 \p ^ ( 71 , 71 )), we construct the final 
state partition function in any sector of staggered magnetization M s G {—A r2 /2, ...,N 2 /2} according to 


z m b [j] = ex P ( H 3* S f) ‘ 


2-rr 




x S2=±i 


= ^ dX H zL exp (^(*A(—l) 2l+Z2 +j z )S 3 - i\M s ) . 

7r ' S3=±' 


( 65 ) 
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Within each of the M s -sectors, we find that there are 

2 TT 


1 


ZmA o] = I dX [2cos(!)] A exp = ( N 2 


2lT 


0 


A 2 


N 2 \ 


IT+M S J (^ + Ms )\(^-M s )\ 


( 66 ) 


equally probable states. In fact, we can again calculate the spin-spin-correlation function in any M s -sector via 

2tt 


(SMm. = 


1 


y,lvIa 27 tZ Ms [0] 

Again, for x = y we trivially obtain j, whereas in all other cases 


dA II E ^S 3 exp(^iA(-l)^S 3 -TAM s ) 


(67) 


2tt 


(sis i v ) Ma = - 


J dX [2 cos m tan 2 «) exp (-iA M.) = 


( 68 ) 


with (p XtV = (— \Y 1+X2+Vl+V2 = ±1, depending on whether x and y are on the same staggered sublattice or not. 
Accordingly, the spin-spin-correlation function in any M s -sector can be written as 


/ c3 n3\ [iV 4 ( 2M s )~]5 Xi y [A~ (2 M s ) 2 ]ip XjV 

W x^y)M B — 47V 2 (1V 2 - 1) 


(69) 


It has been shown [6(| that the probability p(M s ) of being in a certain M s -sector at to is well-approximated by a 
discrete uniform distribution: 


P(M S ) 


2M s ,E + 1 l M *l — Ms, max 

0 \M a \> M s , max 


(70) 


In fact, M Sjlnax = [A4 S L 2 ] = [0.30743A7 2 ] in the infinite volume limit TV —>• oo at T —> 0. Here, [A] denotes 
the ceiling function which returns the smallest integer not less than A'. For finite values of A, however, there are 
deviations so that M Siinax > [0.30743A 2 ], e. g. M s . max cs [0.39A 2 ] = 25 for N = 8, M Sjmax ~ [0.345A 2 ] = 89 for 
N = 16 and Af Siinax ~ [0.335A 2 ] = 192 for A = 24 [66|. As the measurement-driven real-time evolution does not 
change the A/ s -sector of any ensemble member, this distribution p(M s ) does not change with time. Accordingly, by 
averaging over all allowed M s -sectors and evaluating (l53l) we find for the Fourier modes 


A (3) (p = (7r,7r)) = (|S(p = (tt, 7r))| 2 )(t 0 ) - lM Simax (M s , max + 1) , 

A {3 \p ± (7T,7r)) = 


1 


A 2 - 1 


A 4 

— ^(\S(p=( 7 r,n))\ 2 )(t 0 ) 


1 


A 2 - 1 


A 4 1 


- -M. 


c(Ms,n 


1) 


(71a) 

(71b) 


We emphasize that the relative error of this approximation with respect to the analytic value (1M1) is only of the order 
of a few percent and vanishes in the infinite volume limit A — > oo. 


2. Equilibration and attraction time scales 

We have seen that the behavior of the Fourier modes strongly depends on the existence of conserved quantities in 
the measurement process as well as on the momentum value p. In the following, we concentrate on the most interesting 
cases of the measurement processes O ^ and 0^ 3 \ In these cases, we observed the existence of a non-trivial attractor 
A(t) towards which all equilibrating Fourier modes except the slowest one are driven before the final equilibrium value 
A(p) is reached. As a consequence, the approach towards the final equilibrium is slowed down by the slowest Fourier 
mode in the system. Moreover, larger systems are supposed to equilibrate ever slower due to the presence of Fourier 
modes with ever smaller momentum. 

In the following, we further investigate the equilibration dynamics for measurement process O ^. We emphasize 
that an analogous discussion also applies to measurement process O First, we study the final approach towards 
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equilibrium. To this end, suppose that all Fourier modes except the slowest one (note that there are in total four of 
them on the isotropic square lattice) with momentum value 


l«l = £ (72) 

have already fallen onto the attractor. In fact, the time-dependence of the slowest Fourier mode is well-described by 
an exponential decay 


(\S(pi)\ 2 )(t) = A + C( Pl )exp(^--^—^ , (73) 

with A = N 4 /4(N 2 — 1). Having measured the parameters C(pi) < 0 and r(pi) and taking into account the sum rule 
HMD, we can determine the ultimate approach of the attractor towards the final equilibrium value according to 

(\S(p ? Vi)f)(t) ‘^ 2 A + B(p) exp =A- exp > ( 74 ) 

where 1 2 denotes the time at which the second-slowest Fourier mode with momentum P 2 has fallen onto the attractor. 
This clearly shows that the final equilibration rate l/r(p) of all Fourier modes is determined by the behavior of the 
slowest Fourier mode 


t{p) = r(pi) . 


(75) 


In fact, the second-slowest Fourier mode (again, there are in total four modes with the same value of \p 2 \) is then 
well-described by a sum of two exponential functions 

<|SWI 2 >W = ^-^|«p(-^)+C(b)«p(-5^) , (76) 

one of them describing the approach towards the attractor and the other the flow along the attractor. We can again 
numerically determine the coefficient C(p 2 ) < 0 and the relaxation time 1/T(p 2 ). Employing the sum rule (1541) again, 
we can deduce the attractor according to 


(\S(p^p h2 )\ 2 )(t) ~ 3 A- 


4C(pi) 
N 2 - 5 




4 C(p 2 ) 
N 2 - 9 


exp 



(77) 


where £3 < t 2 denotes the time at which the third-slowest Fourier mode with momentum p% has fallen onto the 
attractor. 

In principle, this procedure can be continued for all momentum modes and therefore allows for the iterative deter¬ 
mination of the attractor 


TV 4 ~ f t \ 

Ait) = HN*-l) + ’E CMa v{-Tg>) ’ (78) 

with T(p\) = r(pi). Most important, the temporal dependence of the attractor is determined by a sum of exponential 
functions with different attraction time scales T(p), which describe the approach of the single momentum modes 
towards the attractor. It is clear, however, that this method becomes increasingly involved for larger momenta due 
to the numerical uncertainties in the data. Moreover, for ever larger momenta the order in which the single Fourier 
modes approach the attractor becomes less clear. For the lowest Fourier modes, however, all these issues are under 
control and we can reliably determine the parameters C(p) and T(p). 

In fact, we can simultaneously fit the slowest Fourier modes i.e. the Fourier modes in the vicinity of the (0, 0)- 
mode for and in the vicinity of the (n, 7r)-mode for respectively - upon taking into account the analytically 
predicted attractor behavior. We show the good agreement of these fits with the numerical data in Fig. [5] The 
precise determination of the parameters further allows us to study the momentum dependence of the time scale T(p), 
which determines the rate at which the different Fourier modes are driven towards the attractor. Fig. [ 6 ] displays the 
momentum dependence of the attraction rates for the different measurement processes, along with the best fit of the 
corresponding data. Most notably, for we find in the small momentum regime 

['yT(p)]- t = ci(a\p\) ri , 


( 79 ) 
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<|S(P)| 2 > <|S(p)| 2 > 




FIG. 5: [Color online] (Heisenberg anti-ferromagnet initial state) Time evolution of certain Fourier modes (|<S'(p)| 2 )(t) for the 
measurement processes O (b (left) and O ® (right). The error bars are of the order of the symbol sizes. The solid lines correspond 
to a simultaneous fit of the twelve slowest Fourier modes upon taking into account the analytically predicted attractor behavior. 
The parameters are 4AT = 512, L = 16a, /3J = 5L/2a = 40, ey = 0.05 and we performed 4 • 10 6 Monte Carlo measurements. 


with Ci = 1.17(2) and ri = 2.05(4), cf. also Fig. lb and the corresponding analysis in (5(|. On the other hand, for 
Ol 3 ) we obtain in the large momentum regime 

[jTtp)]- 1 = c 3 (a\p-(ir,n)\) r3 , (80) 

with C 3 = 1.12(2) and = 1.99(4). Obviously, the attraction rate shows for both measurement processes a quadratic 
momentum dependence. Due to this behavior, we regard the approach towards the attractor as a diffusion process of 
the conserved quantity. Larger systems are supposed to equilibrate ever slower due to the presence of Fourier modes 
with ever smaller momentum (1721) . In the infinite volume limit TV —» oo, the momentum variable becomes continuous 
and the equilibration rate d79l) and (l80l) . respectively, becomes arbitrarily small. 

This is contrasted by the momentum dependence of the equilibration rate for the measurement process O ( 2 ) which 
does not conserve any of the Fourier modes. In this case, all Fourier modes equilibrate very quickly as there is no 
additional constraint which could delay the approach towards the new equilibrium state. 


tyT(p)] 1 


lr Tip)]" 1 




a|p-(7r,jr)| 


FIG. 6 : [Color online] (Heisenberg anti-ferromagnet initial state) Momentum dependence of the attraction rate 1/T(p) for the 
measurement processes and O 1 - 2 ' 1 (left) as well as O ( 3 ) and O^ (right). The data points for the measurement processes 0 ( b 
and C)( 3 ), respectively, are calculated by a simultaneous fit of the twelve slowest Fourier modes upon taking into account the 
analytically predicted attractor behavior. On the other hand, we performed independent exponential fits for the measurement 
process 0^ 2 \ The parameters are as in Fig. [5] 
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FIG. 7: [Color online] (Heisenberg anti-ferromagnet initial state) Time-dependent Binder ratio B 4 (t) for the measurement 
process 0 (2 ^ for different values of L/a along with variable /3J = 2>L/Aa. The error bars are of the order of the symbol sizes 
and the lines are included to guide the eye. We emphasize that the results for the measurement process O ^ are essentially 
identical, cf. also Fig. 2b in [561 ]. 


3. Non-equilibrium phase transition 


In the anti-ferromagnetic Heisenberg model, at low temperatures T —> 0 according to (IM1) we find the behavior 


(M s 2 )(f 0 ) _ (|S(p = (?r,7r)| 2 )(f 0 ) 2 

L 2 L 2 


(81) 


indicating spontaneous symmetry breaking of the SU( 2) spin symmetry. During the Lindblad evolution, this order is 
then destroyed by the measurement processes O W and , respectively, whereas it is conserved for the measurement 
process O^. In the first two cases, the system is driven towards the final states (l6lT) and (linJl) . respectively, for which 
(M 2 )/ L 2 becomes volume-independent, indicating that the SU{ 2) spin symmetry is restored again. In the infinite 
volume limit, the system evolves from an ordered state with a finite order parameter density to a disordered state 
where the order parameter density vanishes. Accordingly, the system must undergo a phase transition at some point, 
however, it is not clear a-priori whether this transition occurs at a certain instant of time or rather takes an infinite 
amount of time. In fact, since the Lindblad evolution drives the system far out-of-equilibrium, this phase transition 
is not expected to fall into any of the standard dynamical universality classes [67j. 

In order to study the temporal behavior of the phase transition, we display the real-time evolution of the Binder 
ratio [H, 69|] 


B 4 (t) = 


[(M*)(t )] 2 


(82) 


in Fig. [TJ Typically, for systems in thermal equilibrium plotting B 4 versus the temperature for different system sizes 
produces curves that intersect each other close to the phase transition temperature. Doing the same thing for the 
real-time evolution, we observe that the various finite-volume curves for B 4 (t) do not intersect each other. Moreover, 
we find that their inflection points move to ever later times with increasing volumes. We interpret this result as 
saying that the phase transition does not occur at any specific instant in time but takes rather an infinite amount of 
time before it is completed. It is quite remarkable that the results for the measurement processes O ^ and O^ are 
essentially identical even though their conservation properties and loop-cluster rules are completely different. This 
indicates that the destruction of the anti-ferromagnetic order is rather insensitive to the specific dissipative process. 

In order to further study the dynamics of the phase transition, we investigated the time-dependence of the staggered 
magnetization density A4 s (f). Based on finite-volume chiral perturbation theory (1641) . we first performed a finite size 
analysis at to to determine the shape coefficients c n up to forth order according to 


(M 2 s ){to) = 


M 2 t {to)L< 


3 

■£< 

n =0 


(«M) 






(83) 


As a reminder, the best numerical values for the used low-energy parameters are A4 s (fo) = 0.30743(l)/a 2 and 
£(fo) = 1.459(3)a [HI |66]- We used finite-volume data up to L = 64 a as well as results from [HI to obtain cq = 1, 
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M s (t)/M s (0) f(t)/f(0) 




FIG. 8: [Color online] (Heisenberg anti-ferromagnet initial state) Left: Time-dependent staggered magnetization density 

_A4 s (t)/.M s (0) for the measurement process O 1 ' 2 ' 1 on a logarithmic plot along with an exponential fit of the data. Right: Time- 
dependent length scale £(t)/£(0) for the measurement process O ^ on a linear plot with the lines included to guide the eye. 
In both cases, the error bars are of the order of the symbol sizes. We again emphasize that the results for the measurement 
process o (1) are essentially identical, cf. also Figs. 2c-d in [5§ |. 


ci = 5.7503(6), C 2 = 16.31(2) and C 3 = —84.8(2). The shape coefficients c n , which depend on the geometry of the 
quantum system, are then assumed to be time-independent as the spatial geometry remains unchanged. Moreover, 
we suppose that the chiral-perturbation-theory-inspired formula remains valid in real-time, so that 




fmv Mim 4 

3 ^ ” V L ) - 3 

71=0 V 7 


3 

£< 

71 = 0 




fm 


(84) 


This then allows us to study the decay of the staggered magnetization in the infinite volume limit up to "ft — 0.35 
with an acceptable y; 2 /d.o.f.. The decay rates can then be estimated from an exponential fit 

A4 s (f) = A4 S (0) exp (-t/r) , (85) 

with an inverse decay rate yr = 0.240(2), as shown in Fig. [ 8 ] Again, this suggests that the phase transition is 
completed only after an infinite amount of time rather than after a finite time interval. Additionally, we find that the 
length scale £(!)/£( 0 ) increases with time, which can be attributed to a decrease of the spin stiffness p s . 


IV. CONCLUSIONS & OUTLOOK 

We investigated the real-time evolution of large strongly-coupled 2-dimensional systems of quantum spins s = ^ 
whose dynamics are entirely driven by measurements of neighboring spin pairs. We considered different thermal initial 
states at low temperature corresponding to the anti-ferromagnetic or ferromagnetic Heisenberg model as well as the 
XY-model. The dissipative time evolution, described by the Kossakowski-Lindblad equation, destroys the long-range 
order of the initial states and drives the systems to a new equilibrium with only short-range correlations. To deepen 
the understanding of the real-time dynamics, we studied different measurement processes, which differed from each 
other in their symmetry properties, and derived the corresponding loop-cluster rules. The path integral along the 
Konstantinov-Perel’ contour was then solved with an efficient loop-cluster algorithm. 

We studied in detail the anti-ferromagnetic Heisenberg model initial state for three different measurement processes, 
of which two result in the conservation of the (0,0)-mode (uniform magnetization) or of the (n, 7 r)-mode (staggered 
magnetization), respectively. The conservation of any of the Fourier modes leads to a drastic slowing down of the 
equilibration process due to the occurrence of a non-trivial attractor, towards which all equilibrating Fourier modes 
except the slowest one are driven before the final equilibrium value A(p) is reached. As a consequence, the approach 
towards the final equilibrium is determined by the slowest Fourier mode in the system. In fact, a detailed study of 
the relaxation rate 1/T(jp) at which the different Fourier modes are driven towards the attractor exhibits a quadratic 
momentum dependence. As a consequence, larger systems are supposed to equilibrate ever slower due to the presence 
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of Fourier modes with ever smaller momentum. Due to the quadratic momentum dependence, we identify the approach 
towards the attractor as a diffusion process of the conserved quantities. 

This behavior is in contrast to the case of a dissipative process which does not conserve any of the Fourier modes. In 
this case, all Fourier modes equilibrate very quickly as there is no additional constraint (besides an ever-present sum 
rule) which could delay the approach towards the final equilibrium state. As a consequence, we do neither observe 
slow momentum modes nor a non-trivial attractor towards which the Fourier inodes are driven. Qualitatively, the 
same behavior is also found for initial states corresponding to the ferromagnetic Heisenberg model and the XY-model. 

Studying the Heisenberg model at low temperatures, we find that (M 2 )/L 2 (anti-ferromagnet) or (. M 2 )/L 2 (ferro- 
magnet) are proportional to L 2 , indicating spontaneous symmetry breaking of the underlying SU( 2) spin symmetry 
at zero temperature. Similarly, the large value of ( M 2 )/ L 2 ~ L 2 in the ferromagnetic XY-model indicates the quasi 
long-range order below the Kosterlitz-Thouless transition temperature. After driving the system with various dissi¬ 
pative processes to its new equilibrium state, we find a disordering of the systems corresponding to the restoration 
of the spin symmetry. For the cases under consideration, we find that the staggered magnetization density (anti¬ 
ferromagnetic initial state) or the uniform magnetization density (ferromagnetic initial state) shows an exponential 
decay in time. This suggests that a non-equilibrium phase transition does not occur at any finite instant in time but 
is rather completed only after an infinite amount of time. 

A drawback of the current study is the fact that the real-time dynamics is entirely driven by the measurement 
process whereas the influence of the unitary Hamiltonian evolution has been disregarded in order to enable Monte 
Carlo importance sampling. Nevertheless, its study is still worthwhile as this idealized dynamics may become realizable 
in optical lattice experiments with ultracold atoms in the future. In fact, the main obstacle for including the effect 
of the Hamiltonian so far has been the complex weight problem which arises from the unitarity of the time-evolution 
operator. We are currently exploring whether the meron-cluster idea, which has already been successfully employed 
in solving severe sign problems [70j, [7l|, is also applicable for solving the Lindblad dynamics including a Hamiltonian. 

However, even though we did not explicitly take into account the Hamiltonian in our numerical simulations, we 
can still make statements about the long-time behavior of the full Lindblad evolution, i. e. the combined dissipative 
Markovian and unitary Hamiltonian dynamics. Regarding the Heisenberg model, for instance, the unit density matrix 
corresponding to © is a stable T —> oo fixed point for a large class of measurement processes, including O^ and 
O®. For the measurement process 0^\ however, the system is still supposed to be driven towards the fixed points 
corresponding to GH) or m , respectively. In this respect, the inclusion of the Hamiltonian could only change the 
characteristic time scales of the problem but not the fixed point structure. 
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Appendix A: Loop cluster algorithm 

In order to simulate the real-time dynamics, we employ an efficient multi-cluster algorithm [5.3l.l5 ll| . In this algorithm, 
local stochastic decisions, which are determined by cluster rules, result in non-local changes of worldline configurations. 
Compared to other Monte Carlo procedures, this allows for a drastic reduction of autocorrelation times. Details on 
cluster algorithms can be found in [73 - l74| . 


a. Cluster rules: Euclidean-time branch 


The first model under consideration is the Heisenberg model, whose Hamiltonian is given by 

H = jY j S x -S v . (Al) 

(xy) 

The model is anti-ferromagnetic for J > 0 and ferromagnetic for J < 0. To construct the partition function, we split 
the Hamiltonian into commuting pieces and discretize Euclidean time by performing a Suzuki-Trotter decomposition. 
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TABLE II: Euclidean-time cluster rules for the anti-ferromagnetic and ferromagnetic Heisenberg model. Solid/dashed lines 
denote binding together parallel/anti-parallel spins, respectively, in the same loop-cluster. 


Taking the 3-direction as the quantization axis and denoting the eigenvalues of S% by s x £ {{,}.}, the non-vanishing 
plaquette weights W[s x ,kSy,kiSx,k+ iSy,fc+i] on a bipartite lattice are calculated according to 

W(l) = W[tt,n] = fm,U] = exp(-^) , (A 2 a) 

W(2 ) = W[t4-,t4-] = = exp cosh (^j , (A2b) 

W( 3) = W[U4t] = = exp Sinh , (A2c) 

where e = 1 /N T determines the lattice spacing in the Euclidean time direction. The loop-clusters are then generated by 
assigning bonds in the following denoted by A or A' (vertical parallel or anti-parallel), B or B' (horizontal parallel or 
anti-parallel) and C or C (diagonal parallel or anti-parallel) - with certain probabilities pi j to any spin configuration. 
Here, i £ {1,2, 3,4} corresponds to the spin configuration as introduced in (IA2I) whereas j £ {A, A', B, B' ,C,C'} 
denotes the bond to be chosen. For the anti-ferromagnetic Heisenberg model, the only non-vanishing probabilities are 


Pi,A =P3,B’ = 1 , 
P 2 ,a = 1 — tanh 

P 2 ,B’ = tanh 


eJU\ 

2 ) ’ 


m 


(A3a) 

(A3b) 

(A3c) 


We find that spins which are parallel in the vertical (Euclidean time) direction or anti-parallel in the horizontal 
(spatial) direction are bound together in loop-clusters. On the other hand, the probabilities for the ferromagnetic 
Heisenberg model are given by 


Px,A = exp 



Pi,c = exp 
P2,A = P3,C = 1 ■ 



(A4a) 

(A4b) 

(A4c) 


In this case, only parallel spins are bound together either vertically or diagonally. The cluster rules for the Heisenberg 
model are summarized in Table El 

The second model Hamiltonian under consideration is given by 

H = -JY j (Slsl + SlS 2 y ) =- J -Y. (# S v + S xS+) 

(xy) (xy) 


(A5) 
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corresponding to the (ferromagnetic) quantum XY-model for J > 0. Performing the same steps as for the Heisenberg 
model but taking the 1-direction as the quantization axis and denoting the eigenvalues of Si by s x € {—>,•<—}, on a 
bipartite lattice we obtain the following non-vanishing plaquette weights 


W( 1) = W[-> = W[<- <-1 -] = exp f \ cosh f , 


W(2 ) = W[-t -»■ t— ] = W[< ->■, ->■] = exp 


V 4 J V 4 

, (e/3J 

-—) C0Sh {— 


W(3) = W[—< - >} = W[< -*,—>■«—] = exp sinh , 

W (4) = W[—> —*, <— <— ] = W[<— = exp s i n h 1 


(A6a) 

(A6b) 

(A6c) 

(A6d) 


The corresponding cluster rules, which bind together only parallel spins and which are summarized in Table IIII1 are 
then determined by 


Pi,A = exp ^ 
Pi,b — tanh 


e/3J 

ep£ 

e(3J 


Pi,c = exp ^-— 

P2,A = P3,C = Pl,B 



^ tanh 
= 1 . 



(A7a) 

(A7b) 

(A7c) 

(A7d) 


b. Cluster rules: Real-time branch 

To derive the cluster rules for the real-time branch of the Konstantinov-Perel’ contour, we have to consider m 
for the various measurement processes. More specifically, we have to further investigate the averaged measurement 
results mi), m and (1451) . For all these measurement processes, it follows from (l47l) that the spin configurations and 
loop-clusters that contribute to the path integral are identical on both branches of the real-time path, s x ,k = s' x k . 
Accordingly, we can restrict our attention to deriving the cluster rules on the forward branch of the contour, which 
are then summarized in Table EY1 

For the measurement process corresponding to the total spin O^ 1 ' 1 = S 2 of two neighboring quantum spins, we find 
(s X} k— lSy.fe— \S x k—\Sy k- 1 1 (1 — £7)1 ® 1 + C'yPk \ s x,k s y,k s x ,k s y,k) = 

\S x ,k-lSy,k-l I (1 — £7)1 + e 7 -Pl \S x ,kSy,k) = 2" J ^s Xtk - 1 ,s Xtk S Syk _ 1 ,Sy t k + -^fis x , k - 1 ,Sy, k fisy, k -i,s x , k ■ (A8) 

Employing the same notation as before, this corresponds to cluster rules for binding together parallel spins both 


plaquette 

(1) 

—>• —> 

or 

<- •<— 


(2) 

•<— 

—> 

or 

—» 

(3) 

4- -> 

—>• <- 

or 

—»■ «— 

<- —»■ 

(4) 

<— •<— 

—»■ -> 

or 


XY bonds 

, | | 

B tanh(^) 

C e~~^~ tanh (^y-) 

* 1 1 1 

c X ‘ 

B 1 


TABLE III: Euclidean-time cluster rules for the XY-model with the 1-direction taken as the quantization axis. Solid lines 
denote binding together parallel spins in the same loop-cluster. 
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TABLE IV: Real-time cluster rules for the different measurement processes O ^ with i £ {1,2,3}, with the 3-direction taken 
as the quantization axis. Solid/dashed lines denote binding together parallel/anti-parallel spins, respectively, in the same 
loop-cluster. 


vertically (real time) and diagonally 

Pi,a = 1~ Y , 

(A9a) 


£7 

Pl ’ c = T ’ 

(A9b) 


Pl,A = P3,C = 1 • 

(A9c) 

For the measurement process = S x S y , on 

the other hand, we obtain 


(sx,k-\Sy,k-is' Xtk _ 1 s' y k _ 1 \{l - e7)l 01 + eiPk\s x ,kSy,ks' Xtk s' ytk ) = 


(s x ,k-iSy,k-i (1 - ey)l + eyPy s Xt kS v ,k) = 

(, _ £7\ r , r 

\ 2 / Sx ' k ~ 1 ,Sx ’ k ® S y,k-l i s y,k ' 2 s x,k-l, — S x ,k^S yj k-l, — S y ,k * 

(A10) 


The cluster rules for this process are more intricate. Most notably, this measurement process allows for plaquette 
configurations of the form IT (4) = IT[tti4-4'] = which had a vanishing plaquette weight in the Euclidean 

time branch of the Heisenberg model. Accordingly, the cluster rules which bind together parallel spins both vertically 
and diagonally whereas anti-parallel spins are bound together in the diagonal direction are determined by 


Pi,A = P2,A = 1 0 

(Alla) 

2 — eq 

£7 

Pl,C = P2,C' = 0 r 

(Allb) 

2 — £7 

P3,C = _P4,C" = 1 • 

(Allc) 

Finally, the measurement process O^ = 5+<S'+ + Sy S~ results in 



(sx,fe— lSy t k— l s x,k—l S y,k—l | (1 €q)l ® 1 + C^/Pk \ s x,kSy,kS x ^S y j^ — 

1 ly J ^Sx,k-l,Sx,k^Sy t k~l,Sy t k 4 7y^S Xt k~l,-Sy t k^Sy t k-l,-S X ,k ' 

(A12) 




23 


In fact, we find a vanishing plaquette weight W{ 3) = W[t = W / [4-tit4-] = 0 on the real-time branch for this 

process. This is also reflected in the cluster rules, which bind together parallel spins vertically whereas anti-parallel 
spins are bound together diagonally 


Pi,A = Pa,C' = 1 , (Al3a) 

P2,A = 1 - y , (A13b) 

P 2 ,C = y • (Al3c) 

Concluding, we remark that these cluster rules are dedicated to measuring Ow with i £ {1,2,3} of spin systems 
which have been quantized in the 3-direction. We consider, however, the XY-model for spins which are quantized 
in the 1-direction. Accordingly, these cluster rules then correspond to different measurement processes: In fact, the 
cluster rules (IA9I) still correspond to measuring the total spin S 2 . On the other hand, the cluster rules (IA11I) are then 
related to measuring S 2 S 2 or S^Sy, and the cluster rules (IA13I) correspond to the measurement process S 2 S 2 — S^Sy. 


Appendix B: Heisenberg ferromagnet initial state 


After discussing the measurement-driven real-time dynamics in the anti-ferromagnetic Heisenberg model in great 
detail in Sec. IIII A1 we now consider an initial density matrix po corresponding to the ferromagnetic Heisenberg model 
(ED, where we choose again the 3-direction as the quantization axis. The ensemble of initial states is then prepared by 
means of the Euclidean-time cluster rules (IA4D . We then again study the time-dependence of the Fourier modes (1551) 
for the different measurement processes with z S {1,2,3}. As the measurement processes are independent of the 
chosen initial state, we again find the conservation of the (0, 0)-mode - which now corresponds to the order parameter 
of the system - for O ^ and the conservation of the (7r,7r)-mode for . On the other hand, the measurement 
process O*- 2 -* does again not conserve any of the Fourier modes. 

Applying the same methods as before and assuming that the system is prepared at low temperature T —> 0 at 
the initial time to, we can again calculate the final equilibrium values analytically. The ground state has total spin 
M = N 2 /2 and is (2A4 + l)-fold degenerate, so that we find 


A {1) (P= (0,0)) = {\S{p = (0,0))| 2 }(f o ) = • M( ^ + 1) = jV2( ^ + 2) (Bl) 

for the measurement process . Due to the fact that the (0,0)-mode is exactly conserved during the time evolution, 
the system is driven to a final equilibrium ensemble for which the 3-component of the uniform magnetization is 
uniformly distributed among the 2A4 + 1 sectors with M £ In fact, the corresponding spin-spin- 

correlation function in any sector of M is given by 


( 5 , 


3 Sy)M = 


[A^ 4 - (2 M) 2 ]S x , y - [IV 2 - (2M) 2 
4 N 2 {N 2 - 1) 


(B2) 


Evaluating (1551) and averaging over the 2A4 + 1 sectors of uniform magnetization then gives 


A«(p=(0,0)) 

A«(p^(0,0)) 


N 2 (N 2 + 2) 
12 
N 2 

~G~ ' 


(B3a) 

(B3b) 


As these values also correspond to their equilibrium values at the initial time t 0l we do not observe any interesting 
dynamics upon applying measurement process to the ferromagnetic Heisenberg model at T —> 0. 

For the measurement process 0^ 2 \ on the other hand, the final equilibrium corresponds again to a state where each 
of the 2 n spin states is equally probable. Accordingly, all Fourier modes again take the same final value 


A^(P) = —. (B4) 

We find that all Fourier modes show a rapid equilibration towards these new equilibrium values, as shown in Fig. [9] 
Note again the overshooting of several Fourier modes which can be traced back to (1551) but the absence of a non-trivial 
attractor. 
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FIG. 9: [Color online] (Heisenberg ferromagnet initial state) Time evolution of certain Fourier modes (|S(p)| 2 )(t) for a short 
time interval for the measurement processes O ^ {left) and O ® {right). The error bars are again of the order of the symbol 
sizes and the lines are included to guide the eye. We initialize a ferromagnetic Heisenberg model at low temperatures such 
that the (0, 0)-mode is very large whereas the ( 7 r, 7 r)-mode is small. The remaining parameters are 4 N T = 512, L = 16a, 
/ 3J = 5L/2a = 40 and £7 = 0.05. In order to determine the time evolution, we performed 10 6 Monte Carlo measurements. The 
horizontal lines corresponds to the analytically derived final equilibrium values (IB4I) and (IB6bl) . respectively. 


Finally, for the measurement process O ^ we cannot use the same reasoning as in the anti-ferromagnet as the 
probability p{M s ) of being in a certain M s -sector at to is unknown. The numerical results, however, indicate that the 
conserved (7r, 7r)-mode takes the value 


A [3 \p = (tt, tt)) = (|S(p = (7T,7r))| 2 )(t 0 ) ^ 


N 2 


(B5) 


Moreover, as it is still true that all Fourier modes (p ^ ( 7 r, 7 r)) reach the same final equilibrium value, we conclude 


N 2 

A (3) {p = (7T,7r)) = (I S(p= (7T,7r))| 2 )(t 0 ) ^ -g- , 


H (3) (p ± (7T,7r)) = 


1 


N 2 - 1 


N 4 


- <1 S(p= (7T,7r))| 2 )(to) 


N 2 {3N 2 - 2) 
12{N 2 - 1) ' 


(B6a) 

(B6b) 


Due to the appearance of slow Fourier modes in the vicinity of the ( 7 r, 7r)-mode, we again find a non-trivial attractor 
for the dynamics which is driven by measurement process 0^ 3 \ as shown in Fig. [9] As a consequence, we can again 
simultaneously fit the Fourier modes in the vicinity of the {n, 7r)-mode for measurement process Cb 3 ) upon taking into 
account the analytically predicted attractor behavior. The attraction rates, at which the different Fourier modes are 
driven towards the attractor A{t), are then again found to depend quadratically on the momentum value (15U1) with 
c 3 = 1.22(3) and r 3 = 2.04(5). 

Finally, we studied again the time-dependence of (\S(p = (0, 0))| 2 )/L 2 , whose large value (IBID ~ 0{L 2 ) indicates 
spontaneous symmetry breaking of the SU{2) spin symmetry at zero temperature. On the other hand, after driving 
the system with the measurement processes and to its new equilibrium state, we find that it becomes 
volume-independent suggesting the restoration of the SU{2) spin symmetry. As a matter of fact, the Binder ratio 
B^{t) shows a very similar behavior as for the anti-ferromagnetic initial state, which indicates that the phase transition 
is again completed only after an infinite amount of time. 


Appendix C: XY-model initial state 

In this section, we briefly discuss the dynamics for an initial density matrix po corresponding to the ferromagnetic 
quantum XY-model (IA5I) . We note that the anti-ferromagnetic and ferromagnetic XY-model are unitarily equivalent 
on a bipartite lattice in the absence of a magnetic field. In order to have direct access to the order parameter of the 
system, we now choose the 1-direction as the quantization axis. The ensemble of initial states is then prepared by 
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FIG. 10: [Color online] (XY-model initial state) Time evolution of certain Fourier modes (|5'(p)| 2 )(t) for a long time interval 
for the measurement processes O p) (left) and O ® (right). The error bars are of the order of the line width. We initialize the 
XY-model at low temperatures such that the (0, 0)-mode is large whereas the (-7T, 7r)-mode is small. The remaining parameters 
are 4 N r = 512, L = 16a, /3J = 5L/2a = 40 and ty = 0.05. In order to determine the time evolution, we performed 10 6 Monte 
Carlo measurements. The horizontal lines corresponds to the final equilibrium values (IC4bl) and (IC5bll . respectively. 


means of the Euclidean-time cluster rules (1A71) . Due to the different quantization axis, the Fourier modes are now 
defined according to 


S(p) = '^2 exp (* pa; ) S x = exp + W2X2) si. 

X X 


(Cl) 


We also note that the measurement processes have a slightly different meaning than before as already discussed 
previously. Nevertheless, we have again the conservation of the (0, 0)-mode for O ^ and the conservation of the 
(7r,7r)-mode for O^ 3 ). On the other hand, the measurement process O ^ does again not conserve any of the Fourier 
modes. 

Preparing an initial state at low temperature T —¥ 0, we study the time-dependence of the Fourier modes for the 
different measurement processes. For this initial state, finite-volume chiral perturbation theory predicts 

(\S(p=(0M 2 )(t O )^-^J2 a 4^) > (C2) 

where the first three coefficients are given by a 0 = 1, ai = 0.45157 and a 2 = 0.030793 [fi^ - [64l| . The low-energy 
parameters are the magnetization density A4 = 0.43561(l)/a 2 , the spin stiffness p = 0.26974(5)J and the spin-wave 
velocity c= 1.1347(2 )Ja [zE [zl]. 

Again, the final equilibrium for the measurement process O^ corresponds to a state where each of the 2 N spin 
states are equally probable so that all Fourier modes take the same value 


A^(p) = ^ . 


(C3) 


As in the Heisenberg model, we then find that all Fourier modes show a rapid equilibration towards these new 
equilibrium values without generating a non-trivial attractor. 

For the measurement processes and O^ 3 ), however, we are not able to predict the final equilibrium values A(p) 
analytically. Nevertheless, it is still true that all Fourier modes except the conserved one converge towards the same 
final equilibrium value. Accordingly, we have for the measurement process Cb 1 ' 


A^\p = (0,0)) = (\S(p = (0,0))| 2 )(t 0 ) , 


A {1) (p^( 0 , 0 )) = 


1 


N 2 - 1 


IV 4 


-(|S(p=(0,0))| 2 >(to) 


(C4a) 

(C4b) 
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M(t)/M(0) 




FIG. 11: [Color online] (XY-model initial state) Left: Time-dependent Binder ratio B±(t) for the measurement process 0 (3) 
for different values of L/a along with variable /3J = 3L/Aa. The lines are included to guide the eye. Right: Time-dependent 
magnetization density A4(t)/A4(0) for the measurement process O ® on a logarithmic plot along with an exponential fit of the 
data. In both cases, the error bars are of the order of the symbol sizes. We emphasize that the results for the measurement 
process O (2) are essentially identical. 


with (|S(p = (0,0))| 1 2 )(t o ) being given by (IC2I) . On the other hand, for the measurement process O^ we find 


A (3 \p = (7T,7r)) = (\S(p = (7T,7r))| 2 )(fo) , 


A (3 \p ± (7T,7T)) 


1 fiV 4 5 6 7 


N 2 - 1 


4 


(\S(p= (7T,7r))| 2 )(t 0 ) , 


(C5a) 

(C5b) 


where (\S(p = (7r, vt))| 2 )(to) needs to be determined numerically. In both cases, we again observe slow Fourier modes 
in the vicinity of the conserved Fourier mode, as shown in Fig. 1101 As a consequence, a non-trivial attractor A(t) is 
formed towards which all equilibrating Fourier modes except the slowest one are driven before the final equilibrium 
value A(p) is reached. The attraction rates, at which the different Fourier modes are driven towards the attractor, 
are then again found to depend quadratically on the momentum value (1 79 II with Ci = 1.13(2) and r\ = 2.04(4) as well 
as (l80l) with C3 = 1.19(3) and r 3 = 2.02(5). 

Finally, we investigate the time-dependence of (\S(p = (0, 0))| 2 )/L 2 , whose large value (IC2I) ~ 0(L' 2 ) is a con¬ 
sequence of the quasi long-range order in the two-dimensional XY-model below the Kosterlitz-Thouless transition 
temperature [77 j. We find that it becomes volume-independent after driving the system with the measurement pro¬ 
cesses O ( 2 ) and O*- 3 ) to its new equilibrium state, corresponding to a complete disordering of the spin system. In 
order to study the transition between these two different phases, we again investigated the Binder ratio 


Bi{t) 


(M 4 )(t) 

[<M 2 )(f)P 


(C6) 


and performed an accurate finite size analysis like for the anti-ferromagnetic Heisenberg model. As shown in Fig. 1111 
the various finite-volume curves for B±[t) do not intersect each other but their inflection points move to ever later 
times with increasing volumes. Moreover, the magnetization density shows an exponential decay in time with an 
inverse decay rate yr = 0.242(2). All these observations indicate again that the phase transition does not occur at 
any finite point in time but is rather completed only after an infinite amount of time. 
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